setwd("C:/Users/hjiang/Documents/my project/breastCancerScreen/data")

bcscreen = read.csv("bcscreen.csv", stringsAsFactors = FALSE)

bcscreen$cancer=bcscreen$nurseRefer=bcscreen$radRefer=bcscreen$anyRefer = 0
bcscreen$cancer[bcscreen$CADETAIL<=3] = 1
bcscreen$cancer[bcscreen$intervalca==1]=1
bcscreen$nurseRefer[bcscreen$NRSREFL=="Y"|bcscreen$NRSREFR=="Y"]=1
bcscreen$radRefer[bcscreen$RADREFL=="Y"|bcscreen$RADREFR=="Y"]=1
bcscreen$anyRefer[bcscreen$nurseRefer==1|bcscreen$radRefer==1]=1




## for affiliated center with CBE
##cancer diagnosised by cbe only
#ACBEcancer = AcenterWithCBE[AcenterWithCBE$cancer==1,]
#CBEonly=ACBEcancer[ACBEcancer$nurseRefer==1|ACBEcancer$radRefer==0,]

library(mgcv)
library(MASS)
library(INLA)                                         
library(lme4)
library(glmmBUGS)
library(nlme)
library(R2WinBUGS)

factorBigBase = function(x) {
  factor(x, levels=names(sort(table(x), decreasing=T)))
}


#mammo sensitivity:true-positive: screen-detected cancer
#            false negative: interval cancer+ CBE alone detect

cancer = bcscreen[!is.na(bcscreen$CADETAIL<=3)|!is.na(bcscreen$intervalca==1),]
cancer$response=0
cancer$response[cancer$radRefer==1]=1
cancer$derestrogen[grep("^[[:space:]]*$",cancer$derestrogen)]=NA
res = glmmPQL(response ~ 
                          factor(derestrogen) + 
                          mdegree1 + 
                          odense + 
                          factor(scrnty) + 
                          SCRNAGE +    
                          log(avg_ann_screens) +   ## avg_ann_screenslt4000 + avg_ann_screensge4000 +
                          silength +   ##  silength + 
                          factor(newsiteclass), 
                      data= cancer, random=~1|site/RADID, family= binomial) 
                       
allEstglm = c(fixef(res), "sigma"=sqrt(as.numeric(VarCorr(res)[[1]]))) 
allSEglm = c(sqrt(diag(vcov(res))), NA)

EstCI=cbind("est"=allEstglm, 
                             "sd"=allSEglm, 
                             "LowCI"=allEstglm-1.96* allSEglm, 
                             "UpperCI"=allEstglm+1.96*allSEglm) 


#######################################




forBugs = NULL


#if (nullmodel == T) {
#try(
#forBugs = glmmBUGS(response ~ 1 , effects="site",  family='bernoulli', data=cancer)
#modelFile = "modelNull.bug"
#)
#}else{
#try(


cancer$fderestrogen=factorBigBase(cancer$derestrogen)
cancer$fscrnty=factorBigBase(cancer$scrnty)
cancer$lavg_ann_screens=  log(cancer$avg_ann_screens)
cancer$fnewsiteclass=factorBigBase(cancer$newsiteclass)
forBugs = glmmBUGS(response ~ 
                          fderestrogen + 
                          mdegree1 + 
                          odense + 
                          fscrnty + 
                          SCRNAGE +    
                          lavg_ann_screens +   
                          silength +   
                          fnewsiteclass, 
                      data= cancer , effects=c("site","RADID"), family='bernoulli')
                                                                           
startingVaules = forBugs$startingValues

source("getInits.R")

fromBugs = bugs(forBugs$ragged, getInits, 
  parameters.to.save = names(getInits()),
  model.file="model.bug", n.chain=3, n.iter=3100, 
  n.burnin=100, n.thin=10, 
	working.directory="C:/Users/hjiang/Documents/my project/breastCancerScreen/data", 
	debug=T) 

bugsResult = summaryChain(restoreParams(fromBugs, forBugs$ragged))

bugsResult =  rbind(bugsResult$scalars, bugsResult$betas)

# merge results in


####### for no cancer, reponse is refer or not by rad#############
######## for referreal by radiographer ############
###  specificity= didn't dignosis after normal mamo/didn't dignosed ####
nocancer = bcscreen[bcscreen$cancer==0,]
nocancer$response = 0

nocancer$response[nocancer$anyRefer==0]=1

nocancer$derestrogen[grep("^[[:space:]]*$",cancer$derestrogen)]=NA
nocancerres = glmmPQL(response ~ 
                          factor(derestrogen) + 
                          mdegree1 + 
                          odense + 
                          factor(scrnty) + 
                          SCRNAGE +    
                          log(avg_ann_screens) +   
                          silength +    
                          factor(newsiteclass), 
                      data= nocancer, random=~1|site/RADID, family= binomial) 
                       
allEstglm2 = c(fixef(nocancerres), "sigma"=sqrt(as.numeric(VarCorr(nocancerres)[[1]]))) 
allSEglm2 = c(sqrt(diag(vcov(nocancerres))), NA)

EstCI2=cbind("est"=allEstglm2, 
                             "sd"=allSEglm2, 
                             "LowCI"=allEstglm2-1.96* allSEglm2, 
                             "UpperCI"=allEstglm2+1.96*allSEglm2) 

nocancer$fderestrogen=factorBigBase(nocancer$derestrogen)
nocancer$fscrnty=factorBigBase(nocancer$scrnty)
nocancer$lavg_ann_screens=  log(nocancer$avg_ann_screens)
nocancer$fnewsiteclass=factorBigBase(nocancer$newsiteclass)
forBugs = glmmBUGS(response ~ 
                          fderestrogen + 
                          mdegree1 + 
                          odense + 
                          fscrnty + 
                          SCRNAGE +    
                          lavg_ann_screens +   
                          silength +   
                          fnewsiteclass, 
                      data= nocancer , effects=c("site","RADID"), family='bernoulli')
                                                                           
startingVaules = forBugs$startingValues

source("getInits.R")

fromBugs = bugs(forBugs$ragged, getInits, 
  parameters.to.save = names(getInits()),
  model.file="model.bug", n.chain=3, n.iter=3100, 
  n.burnin=100, n.thin=10, 
	working.directory="C:/Users/hjiang/Documents/my project/breastCancerScreen/data", 
	debug=T) 

bugsResult = summaryChain(restoreParams(fromBugs, forBugs$ragged))

bugsResult =  rbind(bugsResult$scalars, bugsResult$betas)
